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^ I Abstract. A commonly encountered obstacle in indirect searches for galactic dark matter is 

CN ' how to disentangle possible signals from astrophysical backgrounds. Given that such signals 

are most likely subdominant, the search for pronounced spectral features plays a key role for 

indirect detection experiments; monochromatic gamma-ray lines or similar features related 

^f^ ' to internal bremsstrahlung, in particular, provide smoking gun signatures. We perform a 

^D ' dedicated search for the latter in the data taken by the Fermi gamma-ray space telescope 

during its first 43 months. To this end, we use a new adaptive procedure to select optimal 

target regions that takes into account both standard and contracted dark matter profiles. 

The behaviour of our statistical method is tested by a subsampling analysis of the full sky 

p\ . data and found to reproduce the theoretical expectations very well. The limits on the dark 

j^ I matter annihilation cross-section that we derive are stronger than what can be obtained from 

the observation of dwarf galaxies and, at least for the model considered here, collider searches. 

While these limits are still not quite strong enough to probe annihilation rates expected for 

thermally produced dark matter, future prospects to do so are very good. In fact, we already 

find a weak indication, with a significance of S.lu (4.30") when (not) taking into account the 

look-elsewhere effect, for an internal bremsstrahlung-like signal that would correspond to a 

dark matter mass of ~150GeV; the same signal is also well fitted by a gamma-ray line at 

around 130 GeV. Although this would be a fascinating possibility, we caution that a much 

more dedicated analysis and additional data will be necessary to rule out or confirm this 

option. 
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1 Introduction 

While there is little doubt about the very existence of a cold dark matter (DM) component in 
the present composition of the universe, contributing a substantial fraction of il.^ = 0.229 it 
0.015 to its total energy budget [1], the DM nature remains unknown even 80 years after 
Zwicky first postulated a 'missing mass' in the Coma cluster [2]. At the moment, the leading 
hypothesis for a solution to this puzzle are thermally produced, weakly interacting massive 
particles (WIMPs) as expected to appear in many extensions to the standard model of particle 
physics (see, e.g., Refs. [3-5] for reviews or Ref. [6] for a recent pedagogical introduction) - 
though one should keep in mind that the absence of any clear signs for new physics at the 
CERN LHC may render this theoretical prejudice less likely already in a few years from now 
[7]. If DM consists of WIMPs, it can be searched for both at colliders (with missing transverse 
energy as the main signature), in direct detection experiments looking for the scattering of 
WIMPs with the nuclei of terrestrial detectors, or indirectly through the observation of WIMP 
annihilation products in cosmic rays; the obvious advantage of indirect searches (see [8] for 
a recent status review) being that they allow to probe the nature of DM not only locally but 
on cosmological (or at least galactic) scales. 

Gamma rays are for several reasons a particularly suitable detection channel for DM 
annihilation, not the least because they propagate essentially unhindered through the galaxy 
and thus directly point back to the sources. Here, we will focus on another important aspect, 



namely the possibility of sharp spectral features in the annihilation spectra that would allow 
a rather straightforward discrimination from astrophysical backgrounds. Indeed, it has early 
been pointed out [9] that the direct annihilation into photons would lead to the smoking gun 
signature of a gamma-ray line [10-18]. Equally pronounced spectral features near photon 
energies close to the kinematical endpoint at E^ = m^ arise due to internal bremsstrahlung, 
i.e. in the presence of three-body final states containing a photon [19-24], albeit at much 
higher rates because they are not loop-suppressed. Including such features in the analysis, 
rather than following the common approach of looking for featureless spectral templates, can 
significantly increase the sensitivity to DM signals [25, 26] and might even be used to rather 
efficiently discriminate between different DM models (see e.g. Ref. [27]). 

The Large Area Telescope (LAT) [28], the main instrument on board the Fermi Gamma- 
ray space telescope, has an unprecedented sensitivity to gamma rays from 30 MeV to above 
300 GeV. This, together with its large field-of-view, makes it ideally suited for DM searches [29- 
39]; the non-observation of gamma-ray signals from dwarf galaxies, e.g., places the currently 
most stringent bounds on the annihilation rate of WIMPs with masses below around 700 GeV 
[30, 31, 39] (for higher masses, H.E.S.S. observations of the galactic center lead to stronger 
limits [40]). So far, only line-signals have been searched for in the LAT [36, 37] or EGRET [41] 
data; the main purpose of this article is to extend these searches to other spectral endpoint 
features which, as outlined above, are arguably more relevant in most WIMP models. In 
order to look for such templates, we adopt a strategy that is very close in spirit to traditional 
line searches, as recently shown to be very promising [26], and use a new adaptive method 
to identify optimized target regions around the galactic center. As we will show, this leads 
to very competitive constraints on possible annihilation signals. 

For Majorana (but also scalar) WIMPs, the pair of annihilating DM particles approx- 
imately forms a J = state for the small relative velocities v ~ 10~^ expected in the 
galactic halo. In this commonly encountered situation, the annihilation rate into fermionic 
two-body final states // is (usually quite strongly) suppressed by m^/m^ and the next order 
result — with an additional photon in the final state — can be greatly enhanced by a factor 
of ~ (aem/7r)m?^/mp since this only happens for t- and u-channel annihilation mediated by 
the exchange of charged virtual particles, this process is also referred to as virtual internal 
bremsstrahlung (VIB, see Ref. [22] for an extensive discussion). Apart from the expected 
large rates, VIB also results in very pronounced spectral features close to the kinematic end- 
point which resemble a slightly distorted line. Taken together, these two aspects make VIB 
in some sense the most promising DM signature to look for and motivate our choice of mainly 
focussing on related spectral distortions in the astrophysical background. 

In order not to obscure our analysis by the potentially many parameters entering into the 
DM model, we will focus on a simple toy model that corresponds to the minimal extension 
to the standard model where we can expect strong VIB signals. While this model has 
been considered before and in its own right — in particular in connection with electroweak 
bremsstrahlung corrections to the annihilation rate [42-48], but also in other contexts [49, 
50] — let us stress that our analysis can straight-forwardly be applied to any other model 
with similar gamma-ray spectra from DM annihilation; in particular, the main features of 
this model are the same, for our purpose, as for neutralino DM in some phenomenologically 
very relevant regions of the supersymmetric parameter space. 

The structure of this article is as follows. In Section 2, we introduce our toy model 
and comment on how it relates to the more commonly considered case of supersymmetric 
neutralino DM. We describe our method to search for pronounced spectral templates in the 



LAT data in Section 3 and also present our results there. In Section 4, we compare our 
new limits to existing limits from dwarf galaxies, expectations for thermally produced DM, 
collider constraints and limits from cosmic-ray anti-protons. We present our conclusions in 
Section 5. Finally, we provide some additional technical information about our method of 
selecting a target region optimized for the search of DM-related spectral features (Appendix 
A) and how a subsampling analysis of the full sky data can be used as a further test to 
confirm the reliability of our statistical method (Appendix B). 

2 Particle physics scenario 

2.1 Toy model with large virtual internal bremsstrahlung 

We will assume that the DM of the Universe is constituted by Majorana fermions Xj singlets 
under the Standard Model gauge group, which couple to the Standard Model via a Yukawa 
interaction with a scalar rj that is not much heavier than the DM particle. The Lagrangian 
of the model reads: 

^ = -^SM + -^x + -^'7 + -^int • (2-1) 

Here, £sM is the Standard Model Lagrangian. £^ and £^ are the parts of the Lagrangian 
involving only the Majorana fermion x arid the scalar particle r], respectively, and are given 
by 

^x = l^t^^X - -^m^tx , ^2 2) 

£„ = (Z),,77)t(Z)^,?)-myr?, 

where -D^ denotes the covariant derivative. Lastly, £int denotes the interaction terms of the 
new particles with Standard Model fields. 

We will consider in this paper three toy models where the DM particle only couples 
to the right-handed muons, tau leptons or bottom quarks, respectively, via a Yukawa inter- 
action with the scalar t]. We assume the latter to be an SU{2)l singlet in order to avoid 
constraints from electroweak precision measurements. The gauge quantum numbers of the 
intermediate scalar ij are (1, l)i for couplings with the muon or the tau {i.e. t/ is a SU{3)c 
and SU{2)l singlet with hypercharge Y = 1), and (3, 1)1/3 ^^ couplings with the bottom 
quark. Furthermore, to guarantee a coupling to just one generation of fermions we assign ij 
a muon number L^ = —1, a tau number L,- = — 1 or a beauty number B = —1, respectively. 
Then, the relevant part of the interaction Lagrangian reads 

Ant = -yx'^RV + h.c. , (2.3) 

with "^ = fi, T, b. Note that in principle additional couplings of the form W Hifrj and {'t]'r])'^ 
are allowed (where H denotes the Higgs doublet). We will neglect them throughout this work 
since they do not directly influence the gamma-ray signature we are interested in. 

In these scenarios, DM particles can annihilate into two fermions with a velocity- 
weighted annihilation cross-section which can be decomposed into an s-wave and a p-wave 
contribution. The s-wave contribution reads in lowest order of the relative center-of-mass 
velocity v [51, 52] 

y'^Nc "^/ 1 



where ^ = (m^/?Ti^)^ parametrizes the mass sphtting between the DM particle x ^^^ the 
t-channel mediator rj, and the color factor N^ is one for muons and taus and three for bottom 
quarks. The p-wave contribution at lowest order in v is [50] 



/ \p-wave 2 V c -L + /^ /o c\ 

(H2-bod. = %8,^2 (1 + ^)4 - (2-5) 

It is important to note that the s-wave contribution to the velocity-weighted annihilation 
cross-section of Majorana fermions is helicity suppressed, by the mass squared of the daughter 
fermion, whereas the p-wave contribution is suppressed by the velocity squared of the galactic 
DM particles today, typically v ~ 10"'^. Therefore, the 2—7-2 annihilation cross-section is 
fairly small, and higher order corrections could play an important role. 

Indeed, it was shown in Refs. [53, 54] that the associated emission of a vector boson 
lifts the helicity suppression in the s-wave contribution to the annihilation cross-section. This 
process was later dubbed virtual internal bremsstrahlung (VIB), which together with photons 
from final-state radiation (FSR) constitute the full internal bremsstrahlung (IB) spectrum 
[22]; the corresponding Feynman diagrams are shown in Fig. 1. The explicit expression for 
the annihilation cross-section into two massless fermions and one VIB photon is [45, 55]: 
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+ l^i±^ + ^^^-^^-^nf^l (2.6) 

/i + 1 2/i V/" + iyJ 

which is usually non-negligible and can in certain instances, in particular for small values 
of /i, be considerably larger than the 2—7-2 annihilation cross-sections Eqns. (2.4,2.5); Qf 
denotes the electric charge of / and rj in units of |e|. We emphasize that throughout this 
paper "3-body process" refers to the VIB process only, whereas "2-body process" refers to 
the helicity-suppressed tree-level process XX ~^ ff plus the FSR photons. When explicitly 
referring to (o"'y)2_body in the following, we will therefore multiply Eqs. (2.4,2.5) by a factor 
of {1 + fdx (IN^^^/dx) , where [20] 

dx TT X \ m'j I 

Furthermore, the energy spectrum of gamma rays produced in the 2—7-3 process 
has a very peculiar shape that allows for an efficient search for gamma rays from inter- 
nal bremsstrahlung. Namely, the differential three-body cross-section, as function of the VIB 
photon energy x = E/m^, is given by [22] 

da a^mV^NcQ} ^ { 2x 



dx 327r2m2 ^ ^ [ (^ + i)(^ + i _ 22;) {n + l-xf 

(^ + l)(^ + l-2x)^ / /. + 1 ^2_ 



2(/i+l-2;)3 \^^+i_2x 

To illustrate the peculiar features of VIB, the energy spectrum of gamma rays that is 
produced per annihilation in our toy model is shown in Fig. 2 for the three different final 
state fermion flavours; for definiteness we assume m^ = 200 GeV and a relatively small 
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Figure 1. Feynman diagrams of the processes that contribute in leading order to the three-body 
annihilation cross-section and produce internal bremsstrahlung. The first diagram very roughly cor- 
responds to VIB, the second and third to FSR (but note that these contributions can be properly 
defined and separated in a gauge-invariant way [22]). 



mass-splitting of ;U = 1.1. The spectra of secondary photons that stem from the subsequent 
decay or fragmentation of the produced fermions are derived using Pythia 6.4.19 [56]. Note 
that in case of bottom-quark final states we also take into account the production of VIB 
gluons following Refs. [48, 57].^ For two-body annihilation, we cross-checked our results 
with the analytical fits from Ref. [58, 59] and find very good agreement. From Fig. 2 it 
is clear that for small enough mass-splittings the gamma-ray spectrum at high energies is 
completely dominated by VIB photons, which show up as a pronounced peak at energies 
close to the dark matter mass. Secondary photons and FSR only become relevant at lower 
energies, or for larger values of /i. In our spectral analysis of galactic center fluxes presented in 
Section 3, we will entirely concentrate on the spectral VIB feature and neglect the featureless 
secondary photons. We will consider the range 1 < /x < 2, because the VIB feature is most 
important in the nearly degenerate case. In this range, the shape of the VIB spectrum is 
almost independent of fi (it becomes slightly wider for larger //), but its normalization can 
vary rather strongly: for fi = 1.1 (// = 2.0), the rate is already suppressed by a factor of 0.55 
(0.05) with respect to the exactly degenerate n = I case; for large n, the rate scales as oc fi"^ 
(whereas the two-body annihilation rate scales like oc ^~^). For comparison with our main 
results, we will also derive limits from dwarf galaxy observations (see Section 4.1); in this 
case we will take into account both VIB and secondary photons. 

2.2 Connection to the MSSM 

Before continuing, let us briefly mention the connection between our toy model and the much 
more often studied case of supersymmetry. The minimal supersymmetric extension to the 
standard model (MSSM) is extremely well motivated from a particle physics point of view — 
leading, in particular, to a unification of gauge couplings and strongly mitigating fine-tuning 
issues in the Higgs sector — and the stability of the lightest supersymmetric particle (LSP) 
is guaranteed by the conservation of i?-parity; if it is neutral and weakly interacting, the 
LSP thus makes for an ideal DM candidate (for a comprehensive and pedagogical primer to 
supersymmetry and the MSSM see e.g. Ref. [61]). 

In most cases, the lightest neutralino is the LSP, and thus a prime candidate for WIMP 
DM [3]. It is a linear combination of the superpartners of the neutral components of the 



^We use throughout the values Os — 0.118 and Qem ~ 1/128 as evaluated at the mass of the Z boson. For 
DM masses m^ = 40 to 300 GeV this approximation affects the VIB photon cross-section at the few percent 
level, and the gluon VIB cross-section by < 20%. 
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Figure 2. Gamma-ray spectrum [N denotes the number of photons produced per annihilation) as 
predicted by our toy model for different final-state fermions, assuming m^ = 200 GcV and a mass- 
splitting of /x = 1.1. Solid lines show the full contribution from three-body final states, including the 
VIB photons close to a; = 1; dotted lines show contributions from the helicity-suppressed two-body 
final states including FSR (in case of muons, the latter is strongly suppressed and not visible on the 
plotted scales). Branching ratios are calculated according to Eqns. (2.4) and (2.6). In case of bottom- 
quarks, we also include contributions from gluon VIB, xx ~^ bbg, following Ref. [48, 57] (dashed line). 
Note that we convolve the spectra shown here with the Fermi LAT energy dispersion as derived from 
the instrument response functions (about AE ^ 10% at around 100 GeV [60]) before any fits to the 
data are performed. 



C/(l) X SU{2) gauge as well as Higgs fields, 
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(2.9) 



and thus a Majorana fermion just like the DM particle in our toy model. As pointed out 
above, the annihilation into fermion-antiferniion pairs // is therefore helicity suppressed in 
the limit of small velocities; this helicity suppression can be lifted if an additional photon is 
present in the final state and annihilation happens via the t-channel exchange of a charged 
particle. In the case of supersymmetry, this can only be achieved through the corresponding 
left- and right-handed sfermions //, and /r which, in the limit of vanishing mf, couple to the 
neutralino and fermions as 
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VLXflfL + VRXfRfR + h.C. 



(2.10) 



where as usual //j/^ = 2(1 i 75)/. Compared to Eq. 2.3, the sfermions thus play exactly 
the same role as t] and the main difference to our toy model is that i) there are two relevant 
scalars for each fermion final state and that ii) the interaction strength ymn is no longer 
a free parameter but uniquely defined by gauge symmetry, and of course the composition of 
the neutralino (see e.g. Ref. [62]): 



2Qf=Fl g 

VL = -j= — gtcmOy^Nii =F ^^12 ■ 

yn = V2QfgtanewNn , 



(2.11) 
(2.12) 



where g is the usual SU(2) couphng and 0w the Weinberg angle. In Eq. (2.11), the upper 
(lower) signs apply if the third component of the weak isospin for / is given by T^ = +1/2 
(^3 = -1/2). 

While the couplings are fixed, the relative importance of the ffj final state to the total 
annihilation rate in the u — >• limit strongly depends, as we have seen, on the mass difference 
fir = (mr/m^)'^ between neutralino and sfermion. In a phenomenological description of 
the MSSM, i.e. without referring to a specific mechanism for supersymmetry breaking, all 
fi? can be treated as essentially free parameters; choosing one sfermion to be much more 
degenerate in mass with the neutralino than all the others thus effectively, for the purpose 
of our discussion, reduces the general MSSM case to our toy model. ^ 

In fact, such a situation is phenomenologically very relevant already in the most minimal 
supersymmetric setup, the constrained MSSM (CMSSM), and appears both in the f- and 
(for neutralino masses larger than what we are interested in here) in the i-coannihilation 
region — which apart from the focus point, funnel and bulk region are the only regions in the 
CMSSM parameter space where the neutralino acquires the correct thermal relic density (see 
e.g. Ref. [63] for a discussion; note, however, that the bulk region now is already excluded by 
LHC data). Correspondingly, it was found in Ref. [22] that VIB can indeed greatly dominate 
the total photon spectrum at high energies in the MSSM in general, but in particular in 
the coannihilation regions of the CMSSM. If only the stau is degenerate in mass with the 
neutralino (as is the case for example in benchmark model BM2 of that reference), the shape 
of the resulting gamma-ray spectrum at high energies is almost indistinguishable from the 
T~^T~'-f case displayed in Fig. 2; the same holds for the more common situation that all 
leptons are rather degenerate with the neutralino (as in benchmark model BM3; small shape 
differences in the high energy spectra can be attributed to slightly different mass splittings). 
BM3 lies with a neutralino mass of m^ = 233.3 GeV in the reach of Fermi LAT, and we will 
comment on this particular benchmark point in light of our results below. 

3 Fermi LAT search for VIB signatures at the galactic center 

The VIB gamma-ray signal produced in our toy model is sharply peaked at energies close to 
the DM mass, very much like gamma-ray lines that are produced in the two-body annihilation 
into 77. For this reason, our analysis will closely follow the spirit of traditional gamma-ray 
line searches [36, 37, 41], as some of us recently proposed in Ref. [26]. The practical advantage 
of DM signals that are very concentrated in the energy spectrum is that searches for these 
signals can be restricted to relatively small energy windows, being only a factor of a few larger 
than the energy resolution of the instrument. For such small energy ranges, it is reasonable 
(and in fact a posteriori justified by our results) to approximate astrophyscial background 
fluxes by a simple power-law; its normalization and spectral index are obtained directly from 
a fit to the data. Since a detailed understanding of the background sources is not necessary in 
this case, it is possible to choose even complex (but very promising) regions like the Galactic 
center as target in the search for a DM signal, as we will do in our present analysis. 



^ Note, however, that y still cannot be treated as a free parameter. Let us also stress that even if choosing 
one /if so small that the annihilation XX ~^ ffl dominates the gamma-ray spectrum at high energies, this 
does not imply that this process (or xX ~^ ff) ^Iso dominates the gamma-ray spectrum at E^ <C m^ or is 
most important in setting the relic density. The effective equivalence between our toy model and the MSSM 
in the limit /x — >■ 1 thus really only refers to the form of the gamma-ray spectrum at energies at and slightly 
below 771 V 



We will here concentrate on DM masses in the range 40 GeV < m^ < 300 GeV. The 
lower end of the mass range is motivated by the LEP constraint on the mass of charged 
exotic particles, which reads rrir^ > 40 GeV (see Section 4.3); since we are mainly interested 
in the degenerate scenario where m^ w ?7i^, this already excludes much lighter dark matter 
particles. For DM masses above 300 GeV, on the other hand, the spectral feature would be 
outside of the nominal energy range of the Fermi LAT. As discussed above, we concentrate 
here on mass splittings in the range 1 < /x < 2, for which the spectral shape of the VIB signal 
is practically independent of /i. 

3.1 Dark matter signal from the galactic halo 

The gamma-ray flux from DM annihilation in the galactic DM halo is given by a line-of-sight 
integral over the DM density squared. 



dJj _ (av) (IN I 2 

dMi ^^' ~ 8^ dE ii.o.^. ^^ 



(0 = 3rZ^3F/ dspi^ir). (3.1) 

X 



Here, ruy^ is the DM mass, {av) the total DM annihilation cross-section averaged along the 
line of sight, dN/dE the energy spectrum of produced gamma rays, and ^ denotes the angle 
to the Galactic center. The coordinate s > runs along the line of sight, and the distance 
to the Galactic center r is given by r{s, ^) = y^r\)^^scos^p + (ssiri^P^, where r^ = 8.5 kpc 
denotes the distance between Sun and the Galactic center. 

We will consider the following generalized Navarro- Frenk- White (NFW [36, 64]) profile 

p^(r) oc o — , (3.2) 

'^''' ' (r/r,)« (1 + r/r,)3-" ' 

normalized to the fiducial value p^ = 0.4 GeV cm~^ at Sun's position [65, 66] and with a 
scaling radius of r^ = 20 kpc. In case of an inner slope of a = 1 this reproduces the standard 
NFW profile. The possible impact of adiabatic contraction [67-70] can be studied in an 
effective way by allowing for larger inner slopes of the profile. We will concentrate on the range 
1 < a < 1.4, which is still compatible with microlensing and dynamical observations [71] 
(traditional adiabatic contraction following Ref. [67] would give rise to even larger values of 
a, see e.g. Ref. [71]). The modified isothermal and Einasto profiles [72-74] are expected to 
give comparable results to the standard NFW profile in searches for line-like features [37]. 

3.2 Event and target region selection 

The gamma-ray data measured by the Fermi LAT is publicly available [75]. The events 
that enter our main analysis are taken from the P7CLEAN_V6 event class. We consider 43 
months of data (from 4 Aug 2008 to 6 Feb 2012), and select front- and back-converted events 
with energies in the range 1-300 GeV. We apply a zenith angle cut of 9 < 100° in order 
to avoid contamination with photons from the earth albedo, as well as the quality cut filter 
DATA_QUAL==1 (all event selection is done using the 06/10/2011 version of ScienceTools 
v9r23pl. 

The target region for our spectral analysis is chosen by using a data-driven adaptive 
procedure with the aim of maximizing the expected signal-to-noise ratio. We stress that 
such an approach is extremely important when looking for spectral features at the statistical 
limit of the detector, and that an inefficiently chosen target region can easily wash out or 
hide a potential signal. Our choice of the optimal target region depends on the adopted DM 
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Figure 3. Target regions used in our spectral analysis (solid black lines). From top left to bottom 
right, Regl, Reg2, Reg3 and Reg4 are respectively optimized for DM profiles with inner slopes of 
a = (1.0, l.f, 1.2, 1.4) as described in the text and in Appendix A. The optimization maximizes 
the signal-to- noise ratio. For comparison, the colors show the expected signal-to- &ac%roMn(i ratio, 
normalized in each case to 1 for the central pixel. 



profile; in order to select it, we estimate the expected spatial distribution of background noise 
in our search for spectral features above 40 GeV by considering the actually measured events 
below 40 GeV. The spatial distribution of signal photons, on the other hand, just follows 
from Eq. (3.1). All details of the method are given in Appendix A. 

We adopt four reference values for the inner slope of the DM profile, a = 1.0, 1.1, 1.2 
and 1.4, for which we obtain the target regions that are shown in Fig. 3 as solid black lines. In 
this plot, the colors encode the expected signal-to-background ratio in different regions of the 
sky, normalized to one for the pixel where this ratio is maximal (note that the actual value 
of this quantity is a factor of 1.9 (3.9, 31) larger for Reg2 (Reg3, Reg4) than for Regl). In 
case of a standard NFW profile with a = 1.0, the target region includes besides the galactic 
center also regions at higher and lower latitudes up to |6| < 70°; for steeper profiles the 
optimal target regions shrink drastically to regions closer to the galactic center. The galactic 
disc is strongly disfavoured in all cases. Southern regions are somewhat preferred, since the 
diffuse gamma-ray emission from our galaxy is not perfectly north/south symmetric. From 
these four regions we extract the measured spatially integrated gamma-ray energy spectrum 
for our subsequent analysis. 



3.3 Spectral analysis 

The search for VIB signatures is done by using the sliding energy window technique discussed 
e.g. in Refs. [26, 36, 37, 41]: we consider for each DM mass m^ in the range 40 GeV < 
m^ < 300 GeV a small energy window that is approximately centered on m^, and hence 
on the position of the expected VIB feature. More precisely, we use the energy range E = 
m^^e~^''^ . . . min[m^e'^'^,300GeV], where the size of the window e varies between e ~ 1.8 for 
m^ = 40 GeV and e ~ 7 for niy. = 300 GeV. The size of the window is identical to the 
values used in Ref. [37] , where it was found to lead to reasonable background fits in context 
of gamma-ray line searches. The position of the window is not exactly symmetric around m^^, 
but slightly shifted towards lower energies as it was suggested for VIB features in Ref. [26] 
in order to increase the sensitivity. We emphasize again that secondary photons, as they 
come from the decay or fragmentation of the fermions, are neglected in our spectral analysis, 
because they become relevant only outside of the energy window that we consider here (at 
least for our toy model). 

For each given mass m^, and within the adopted small energy window, we now fit the 
gamma-ray spectra measured in the different target regions of Fig. 3 with a simple three- 
parameter model: The astrophysical background fluxes are approximated by a power law 
with a free spectral index (1) and normalization (2); the DM VIB signal has only a free 
normalization (3), whereas its mass and the mass-splitting (which we set to /i = 1.1 in most 
of the analysis) remains fixed during the fit. For physical reasons we require the normalization 
of the VIB signal to be positive. 

Technically, we perform a binned analysis of the gamma-ray spectrum measured in the 
different target regions. To this end, we distribute the corresponding measured gamma-ray 
events in a very large number of energy bins (200 per decade). Since the size of the adopted 
energy bins is much smaller than the energy resolution of the Fermi LAT, our analysis is 
essentially identical to an unbinned analysis of the energy spectrum. Each energy bin j then 
contains a number Cj of events. The number of expected events fj,j are obtained by convolving 
our above three-parameter model with the energy dispersion and the exposure of the LAT; 
the resulting fij can then be fitted to the observed counts Cj by maximizing the likelihood 
function C = IljP{cj\fij) with respect to the three model parameters. Here, P(c|/u) is the 
Poisson probability to observe c events when // are expected. Note that the functional form 
of the energy dispersion is directly inferred from the IRF of the P7CLEAN_V6 event class and 
correctly averaged over impact angles and front- and back-converted events, using our own 
software. Exposure maps are derived using the ScienceTools v9r23pl. 

Limits on or the significance of a DM VIB contribution can then be derived by using 
the profile likelihood method [76]. A one-sided 95% C.L. upper limit on the annihilation 
cross-section is obtained by increasing the DM signal normalization from its best-fit value 
until —2lnC increased by a value of 2.71 (while refitting the background parameters). The 
significance of a signal, on the other hand, is derived from the test statistics (TS) 

r5 = -21n /""^^ , (3.3) 

-^best-fit 

where /^best-fit is the likelihood of the best-fit model, and i2nuii the likelihood of the null 
hypothesis (the absence of a DM signal; the null model has only two free parameters). In 
absence of a signal, one expects that the TS follows some x^ distribution. More precisely, 
since the normalization of the DM signal is bounded to be positive, the TS should follow a 
0-5Xfc=o + 0-^Xfc=i distribution (see e.g. Ref. [37]), where x'k=o ^^'^ Xk=i have zero and one 
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Figure 4. Our results for the 95% CL limits on the three-body (VIB) annihilation cross-section of 
XX ~^ f fl: for different values of the DM profile inner slope a. The limits are obtained by a spectral 
analysis of the gamma-ray flux in the corresponding target regions shown in Fig. 3. We assume a 
mass splitting of /i = 1.1. Note that the limits do not directly depend on the nature of the final 
state fermion /, as they are derived from the VIB feature only. The gray cross shows the CMSSM 
benchmark point BM3 from Ref. [22]. 

degree of freedom, respectively.^ This theoretical expectation is indeed very well confirmed 
by a subsampling analysis of the data as we discuss in Appendix B. Hence, if the test statistics 
is measured to be TS for a certain DM mass in a single trial, this would correspond to slightly 
more than \/TSa significance. 

However, since in the present analysis we effectively perform many statistically indepen- 
dent trials when scanning through m^ and analyzing different target regions, the probability 
to find just by chance a statistical fluctuation that mimics a signal is increased; this is known 
as the look-elsewhere effect (LEE). In our case, we approximate the distribution of maximal 
TS values from which the significance of our signature is calculated by 4 x 4 = 16 trials over 
a Xk=2 distribution. Four trials over a Xk=2 distribution come from the scan over tti^ (see 
Appendix B or, for a general discussion, Ref. [77]); the remaining four trials are associated 
with the four target regions. In practice, the significance of the observed signature is then 



for a, where Pixl < x) denotes the 



found by solving P(xi=2 < r5)#*"^'^ = P(xLi < ^^ 

probability to observe a value smaller than x when drawing from a xi distribution. 

3.4 Results 

In our analysis, no VIB signal with a significance of at least 5a were found. Instead, we show 
in Fig. 4 the upper limits at 95% CL on the VIB three-body cross-section, for the different 
values of the inner slope a that correspond to the target regions in Fig. 3. We assumed 
a mass-splitting of /.t = 1.1 for definiteness; limits for other values of /i are very similar 
and will be discussed below (see e.g. Fig. 6). As can be seen from the plot, our limits are 
always stronger than the 'thermal cross-section' of 3 x 10~^^ cm'^ s~^ that is often quoted 
for comparison; in the case of contracted profiles with a = 1.4 they can even reach down 
to values of 10"^^ cm^ s~^ for DM masses ttt-v < 100 GeV. As we will discuss below in 



^The probability distribution function of Xk=o is just 5{TS). For discussions about tlie coverage of confi- 
dence intervals on bounded parameters see Ref. [76]. 
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Figure 5. Left panel: statistical significance for VIB-signal in terms of the TS value, as function of 
rriy^ and for the different target regions shown in Fig. 3. Right panel: fits to data in Reg2 for the 
best signal candidate at to^ = 149 GeV. We show the background-only fit without DM signal as red 
bars. The green bars show the background plus DM signal fit, the blue line the corresponding VIB 
signal flux. In the right panel, we rebinned the data into (9 times) fewer bins than actually used in 
our statistical analysis in order to improve the optical appearance of the figure. Note that the shown 
fluxes are already integrated over the individual energy bins and properly convolved with the LAT 
IRF. 

Section 4.1, our limits are much stronger than what can be obtained from e.g. dwarf galaxy 
observations. For comparison, the gray cross in Fig. 4 shows the CMSSM benchmark point 
BM3 [22], which lies in the coannihilation region and was already discussed above. This 
benchmark point still remains unconstrained by our limits; its rather small cross-section is 
closely related to the requirement that the neutralino is a thermal relic, as we will discuss in 
Section 4.2 below. 

In the left panel of Fig. 5, we show the significance for a VIB-like spectrum as function of 
m^, assuming that /i = 1.1. The different lines correspond to the different target regions. The 
significance is shown in terms of the TS value that was discussed above. We find a possible 
signal candidate at a DM mass of m^ ~ 150 GeV. The indication for a signal is largest 
for the target region Reg2, which corresponds to a = 1.1, and has a nominal significance 
of VTS = 4.3cr. Taking into account the LEE as discussed above, the significance is S.lcr. 
The corresponding fit to the data is shown in the right panel of Fig. 5; the spectral feature 
in the measured flux can be easily recognized by eye. A similar preference for a signal, 
although with less significance, appears also in the other regions Regl, Reg3 and Reg4 (note 
that the fluctuations around 50 GeV are completely within the statistical expectations). TS 
values of zero indicate that for these values of m^j- the data would be best fitted with an 
unphysical negative signal normalization; in this case, the likelihood of the model with DM 
contribution becomes identical to that of the null model because we enforced a non-negative 
signal normalization in our fits. 

We have performed several tests to exclude the tempting DM interpretation of this 
signature, none of which has succeeded so far: By masking out different halves of the signal 
region of Reg2, for example, we find that the signal independently appears in the north, 
south, east and west parts of Reg2 (though with a large scatter in the significances), as 
expected from a DM signal. When shifting the target region away from its position by 
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about 10-20°, on the other hand, the signal disappears completely. This makes a purely 
instrumental effect, which would likely also appear in other sky regions, less likely. To test 
our statistical method, we performed a subsampling analysis of the Fermi LAT data in the 
galactic anticenter region in order to obtain an empirical understanding of the statistical 
behaviour of our TS (see Appendix B for details). We find that — in absence of a signal — our 
TS follows very well a x^ distribution, as theoretically expected. The signature appears also 
in the P7ULTRACLEAN_V6 and P7SOURCE_V6 event classes, and when considering 
back- and front-converted events separately; furthermore, as expected for a real signal, its 
significance has grown with time (i.e. we find smaller TS values when considering instead 
older and older data sets, though with some scatter around the linear trend). The uncertainty 
in the effective area as relevant for searches for line-like features is about 2% [78], which is 
much smaller than the fractional contribution of the observed signature to the flux. 

Despite these encouraging facts, we call for caution when interpreting this signature as 
due to DM annihilation. First of all, with a significance of 3.1cr (when taking the LEE into 
account), it could still simply be an upward fluctuation at the right place. Alternatively, 
the observed signal could be real but due to a yet unidentified astrophysical process, like 
e.g. the inverse Compton scattering of an extremely hard electron component on stellar light 
(see e.g. Ref. [79] for a related discussion of ICS in the context of the Fermi bubbles) — 
though in general, as already stressed several times, it is quite difficult to explain this kind 
of spectral features with astrophysical processes. In this context, it is also worth to mention 
again that our analysis crucially depends on the commonly adopted assumption that the 
background locally follows a power-law, i.e. within each energy window that we consider. 
In principle, it might thus make for an interesting follow-up study to perform a signal fit 
on a more complicated background which contains, e.g., a break in the spectral index that 
could be confused with a signal. We do not expect, however, that such an analysis leads 
to qualitatively different results than presented here because the data itself tells us that 
the single power-law assumption works very well (see Appendix B), and because of the 
sharpness of the observed signature. Even more importantly, we always find a spectral index 
for the background contribution that is consistent with —2.6; this value is expected for the 
production of gamma rays from the collision of cosmic rays with the interstellar medium and 
thus extremely well motivated from astrophysics (see e.g. Ref. [80]). Lastly, our analysis relies 
entirely on the publicly available data, which makes it impossible to take into account all 
known instrumental effects. However, we strongly believe that our actual statistical analysis 
is sound, and not significantly affected by the obvious systematics of the LAT. 

Finally, we note that the best-fit three-body cross-section for a VIB signal in Reg2 is 
(av) = (6.2 lb 1.5 lii'^) X 10"^'^ cm^ s~^ (assuming a = 1.1), the best-fit mass is in the range 
rriy^ = 149 zt 4 j^j^^ GeV; the errors are respectively statistical and systematical.^ This is 
quite a bit larger than what is actually expected for VIB from thermally produced DM (see 
Section 4.2 for a discussion), making a straightforward interpretation of the signal in terms 
of a VIB signal somewhat less appealing. However, let us point out that the signature can 
also well be fitted by a pure gamma-ray line as produced in xx ~^ 77? with a dark matter 
mass of around m,-^ ~ 130 GeV and an annihilation cross-section (av) ~ 10~^^ cm^ s~^ (see 
also Ref. [81] for a dedicated analysis). Finally, we note that the more commonly adopted 
annihilation spectra from e.g. xX ~^ l^^ l^~ or xX ~^ ^^ annihilation are much too flat to fit 
the data and thus cannot be used to explain the signal. 



* Systematical errors stem from uncertainties in the overall effective area (10%) and in the energy calibration 
{t\o%), see Ref. [78]. 
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Figure 6. Comparison of different upper limits on the three-body annihilation cross-section of xx ^ 
fi'^(i~"f, for three reference values of the mass splitting /i. Black lines show 95% CL upper limits that 
come from our spectral analysis of the Galactic center fluxes as shown in Fig. 4, assuming a standard 
NFW profile. Green lines (partially overlapping) show the corresponding limits derived from dwarf 
galaxy observations, taking into account both two- and three-body annihilation channels. Blue lines 
show upper limits on the annihilation cross-section that are derived from requiring that the relic 
density predicted by our toy model does not undershoot the observed value. 
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Figure 7. Same as Fig. 6, but for annihilation into t^t 7. 

4 Comparison with other constraints 

The limits on the cross-section (iTi')3_body as shown in Fig. 4 were derived from VIB photons 
only, according to the spectrum given in Eq. (2.8). We will now discuss these constraints in 
light of other complementary probes, namely limits from gamma-ray observations of dwarf 
galaxies, hmits that can be derived from the thermal production of DM and colhder hmits; 
we finally stress that small values of /i, for the annihilation into 667 final states, can also very 
efficiently be probed by both cosmic ray antiprotons and direct DM detection experiments. 
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Figure 8. Same as Fig. 6, but for annihilation into 667. Note that we include bbg final states when 
calculating the dwarf limits [48, 57]. 

4.1 Limits from dwarf galaxies 

Important limits on the DM annihilation cross-section derive from gamma-ray observations 
of dwarf galaxies [30, 31, 39]. To compare our above results in Fig. 4 with the constraints 
that typically can be obtained from these targets, we will make use of the flux limit presented 
in Ref. [39].^ It reads, in terms of the total annihilation cross-section, {av) < Sir mi/N}^°^ ■ 



5.0 X 10^3° cm^ s 



"x' 1 



Here, N}^°^ denotes the number of photons that are produced 



1 GeV" 

per annihilation in the energy range 1-100 GeV. When calculating N}^"^, we fully take into 
account all photons from three-body as well as the common two-body final states according 
to Eqs. (2.4) and (2.6);^ the relative importance of two- and three-body final states depends 
on fi and the mass of the final state fermions. 

We plot the resulting limits in Figs. 6, 7 and 8 as green lines for difl'erent values of the 
mass splitting fi and for the three diff^erent final state fermion fiavours; in order to allow a 
simple comparison with our result in Fig. 4, we choose to present them in terms of the three- 
body annihilation cross-section (by rescaling them by a factor of ((7t>)3_body/[(<7v)2-body + 
(c^)3-body])- In niost cases the limits depend relatively strongly on the mass splitting fi 
(Figs. 7 and 8); only for annihilation into /i"*"^" the limits remain practically independent of 
fi, since the fiuxes are always dominated by VIB photons (c/. Figs. 2 and 6). In contrast to 
that, our limits from the spectral search for VIB features in the galactic center, as shown by 
the black lines, exhibit only a very weak dependence on the mass splitting parameter /x for 
values of // ~ 1.1-2, since the spectral shape of the VIB radiation changes only mildly in this 
range. 

As long as the mass splitting fi remains small, and for the leptonic final states, our 
spectral search in the galactic center fiuxes for VIB features leads to much stronger constraints 
on our toy model than dwarf galaxy observations, see Figs. 6 and 7. For a large enough mass 



^ Refs. [30, 31] present limits on different common annihilation channels that are however difficult to 
translate into limits on VIB. 

^For simplicity, we here approximated the spectrum of secondary photons from the VIB process with the 
secondary photon spectrum of the corresponding two-body process. This makes the dwarf limits slightly too 
strong at high DM masses (by less than a factor of two). 
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splitting, or in case of colored final state fermions like bottom-quarks, however, the 3-body 
annihilation into VIB photons becomes less relevant and the dwarf limits start to overtake 
our VIB limits (c/. Fig. 8). Note that for colored final states the VIB of gluons becomes 
important [48, 57] and actually dominates the dwarf limits (see also Section 4.4). 

4.2 Thermal production 

In our toy model, self-annihilation in the early Universe is usually dominated by the p-wave 
process XX ~^ f f ■, since the averaged velocity is v ~ C(l), see Eq. (2.5). Any embedding of 
this model in a full theory (like the MSSM, see Section 2.2) would likely increase the total 
annihilation rate. A loose but reasonable upper limit on the coupling constant y can hence be 
derived by requiring that the relic density which follows from the t-channel process xX ~^ ff 
alone does not already undershoot the measured DM density. From this, we can estimate an 
upper limit on the partial annihilation cross-section for xX ~^ ffl ^-s follows: At the freeze- 
out temperature Tf ~ m^/20, the p-wave annihilation channel xX ~^ ff is mildly suppressed 
by a velocity factor of {v'^) = 6Tf/m-^ ~ 0.3; the corresponding annihilation cross-section 
should be equal to or smaller than the thermal cross-section {av)^^^ ~ 3 x 10^^^ cm^ s~^, in 
order to not undershoot the observed DM relic density. However, the s-wave VIB process 
XX ~^ ffl is suppressed by a factor ~ acm/Ti", independently of the temperature T or the 
velocity v. Taking these pieces of information together, one finds an upper limit of roughly 
f3^(f^)th ~ 2 X 10~^^ cm'^ s~^ for the three-body annihilation cross-section (o"f)3_body 
Hence, if we want our toy model to be compatible with the observed relic density (assuming 
thermal DM production and a standard thermal history of the universe), this suggests an 
annihilation cross-section (erf )3_body that is about two orders of magnitude smaller then the 
common thermal cross-section. Note, however, that e.g. destructive interference with s- 
channel diagrams, which do not appear in our toy model, could spoil this argument and thus 
in principle allow for larger rates. 

A full calculation confirms our order-of-magnitude estimates: Neglecting the coannihi- 
lation of x and r], the relic density due to two-body annihilation xX ~^ ff is given by (see 
e.g. Ref. [46]) 

''^'^ -^-'VA y J VlOOGevJ l + /i2 ■ ^^■'> 

The neglect of coannihilation is valid as long as /z > 1.2 [82]; for the value /i = 1.1 and 
m^ = 100 GeV (200 GeV), we have checked with MicrOMEGAs [83] that our Eq. (4.1) 
gives a result just a factor of 2.5 (1.5) larger than the exact value. Requiring now that 
O^/i^ > 0.1, and using Eq. (4.1), we obtain upper limits on the cross-section as shown in 
Figs. 6-8 by the blue lines for different values of /i. These limits largely agree with the rough 
estimates given in the preceding paragraph as long as the mass-splitting is reasonably small. 
For larger mass splittings (or for very small splittings in the regime where coannihilation 
would further reduce the relic density) the limits become stronger. 

From Figs. 6-8 it is apparent that in case of a regular NEW profile with (a = 1) our 
above galactic center limits are still around one order of magnitude away from values of 
the three-body cross-section that are consistent with a thermal relic. However, in case of 
compressed profiles with a > 1, our limits become stronger, as shown in Fig. 4, and for inner 
slopes of the DM halo a ~ 1.2-1.4, they can become sensitive enough to probe the required 
cross-sections of ~ 10~^^ cm'^ s~^ and below. 
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Further enhancements of the signal, without violating our relic density constraints, can 
appear due to the gravitational clustering of DM as predicted by A^-body simulations [74, 
84]. Moreover, the effect of coannihilation could actually be inversed and increase the relic 
density in some cases [85], which again would allow larger annihilation rates while still being 
compatible with the relic density. Finally, we note that the non-thermal production of dark 
matter or a non-standard freeze-out history would invalidate our relic density limits and 
admit larger annihilation cross-sections today (see e.g. Refs. [86-88]). 

4.3 Limits on light charged or colored scalar particles 

Scenarios with light scalars that couple to the Standard Model quarks and leptons are con- 
strained by direct searches at colliders as well as by indirect searches through their effect on 
the quantum corrections to the electroweak precision measurements. 

Collider searches. The most severe constraint on scenarios with light charged particles 
follows from the non-observation of an excess in the invisible decay width of the Z boson at 
LEP, AFinv < 2.0 MeV [89], which excludes the existence of exotic charged scalar particles 
with mass below 40 GeV [90]. Since we are mostly interested in the degenerate case where 
m^ ~ rrijj, we only considered DM masses with m^ > 40 GeV in this work. 

Additional constraints on the mass of charged scalar particles can be derived from 
the non-observation at LEP of an excess over the Standard Model expectations of dilepton 
events with missing energy, generated by the production of charged scalar particles that 
decay into a lepton and an invisible particle (in our framework, the DM particle). The search 
strategy relies on the identification of at least one lepton event as well as on satisfying the 
requirements on isolation and transverse momentum. Hence, the efficiency of the search 
decreases drastically when the DM particle and the intermediate scalar are more and more 
degenerate in mass, which is exactly the region we are interested in. The best constraints on 
such a scenario were published by the OPAL collaboration using 680 pb^^ of e^e~ collisions 
at center-of-mass energies between 192 GeV and 209 GeV; they read [91] 

• "mri ^ 94.0 GeV at the 95% CL for couplings to muons, assuming tti^ — m^ > 4 GeV 
and 100% branching ratio for r] — )• fix 

• 'iT-ry > 89.8 GeV at the 95% CL for coupling to taus, assuming ?n,^ — m^ > 8 GeV and 
100% branching ratio for ij — )■ tx- 

For /_i ~ 1.1, these limits are always satisfied by our toy model. However, much larger mass 
splittings are already partially ruled if the dark matter mass is below ~ 90 GeV. 

Searches for light colored scalar particles have been undertaken both at LEP, Tevatron 
and at the LHC. The non-observation of an excess over the Standard Model background of 
dijet events with missing energy translates into the following 95% CL limits on the mass of 
colored scalars, assuming 100% branching ratio for i] — )■ bx (see also Fig. 4 in Ref. [57] for an 
updated summary): 

• m^n > 875 GeV, assuming mrj-m^> 130 GeV (ATLAS [92]) 

• m^n > 76 GeV, assuming m^^i - m^ > 7 GeV (DELPHI [93]) 

• mrj > 89 GeV, assuming m^^i -m^>8 GeV (ALEPH [94]). 
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Again, as long as // ~ 1.1, our toy model is not affected by these limits; only much larger 
mass splittings are already partially ruled out. 

Oblique parameters. The new scalar particles carry charges under the electroweak group 
and can in principle modify the oblique parameters S, T and U [95, 96]. However, in the 
scenario where the scalar particle is an SU(2)i singlet, none of the three oblique parameters 
receives any exotic contribution. 

Muon g — 2. The interaction which leads to the VIB process XX ~^ /^"""/^"T leads, upon 
closing the DM fermion loop, to a contribution to the muon g — 2. In the toy model presented 
in this work, the contribution from DM particles to the muon g — 2 reads (see e.g. [97]) 

where n = {mrj/m^)'^ denotes the usual mass splitting, which enters the function 

^, , 2 + 3a; - 6a;^ + a;^ + 6a; log X , ^ ^. 

^(^) = 6(^31)^ ■ ^'-'^ 

This contribution should be compared to the 3.2a deviation between the Standard Model 
prediction and the experimental measurement [98], 

<"" - af^ = (29 ± 9) X 10"^° , (4.4) 

which is of opposite sign. Therefore, in our toy model the discrepancy between the theoretical 
prediction of the muon g—2 with the experimental measurement is larger than in the Standard 
Model. If we interpreted the g — 2 anomaly as a statistical fluctuation, the total theoretical 
prediction should still not deviate more than 5a from the experimental value, which implies in 
the limit ;U — )• 1 an upper bound on the coupling of y < 1.7(?7i^/100 GeV). The corresponding 
upper limit on the three-body annihilation cross-section reads 

(^t^)3-body < 2.5 X 10-26 cm3 s~'(^j^^y . (4.5) 

Comparing this to our VIB limits shown in Fig. 4, we find that — within our toy-model and 
the above assumptions — the g — 2 constraints are typically weaker than what we get from 
the spectral analysis of gamma-ray fluxes. Alternatively, additional particles could exist 
that generate a positive contribution to the muon g — 2 which compensates the negative 
contribution from the DM particles, thereby bringing the theoretical prediction closer to the 
experimental value. 

4.4 Anti-proton observations and direct WIMP searches 

In our toy model, the presence of a light colored scalar rj opens up the gluon-VIB channel 
XX ~^ ffg] ill general, this process has a ~ 100 times larger cross-section than xX ~^ f fli 
because instead of aem the strong coupling Ug enters the diagrams [48, 57, 99, 100]. Hence, 
the limits on xX ~^ bbj that can be obtained by constraining the corresponding process 
XX ~^ bbg with cosmic-ray anti-proton observations are quite significant for small values of 
//; depending on the details of cosmic-ray propagation, they can be comparable to or even 
stronger than our gamma-ray limits from the VIB search [48, 57]. Furthermore, xX ~^ bbg 
contributes to the gamma-ray dwarf limits in Fig. 8 [57]. 
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The leptophilic models discussed in this paper also give rise to an antiproton flux through 
the annihilation channel XX ~^ f f^ which, when kinematically open, has a cross-section 
comparable to xX ~^ ffl [42-48]. For these models, our gamma-ray limits are one order of 
magnitude stronger than the antiproton limits in the mass range m^ ~ 100-300 GeV. 

Lastly, for our bottom quark scenario, limits from direct DM searches may become very 
relevant, as the scattering between x ^-^d the detector nuclei via an intermediate rj becomes 
resonant for small values of fj, [101]. 

5 Conclusions 

One of the main challenges that are encountered when searching for signatures from DM 
annihilation in the cosmic gamma-ray flux is the discrimination between a possible signal 
and the astrophysical background. Very pronounced spectral features, like gamma-ray lines 
or photons from internal bremsstrahlung, play a key role in this context because it would be 
difficult to attribute them to astrophysical processes. Previous searches mostly concentrated 
on gamma-ray lines; here, we instead analysed the effects of internal bremsstrahlung which 
is intrinsically more promising due to the larger rates that are expected. 

To this end, we defined a simple toy model (which can be considered as a subset of 
the MSSM particle content) with the important feature to generate an intense gamma-ray 
signature from virtual internal bremsstrahlung. We performed a dedicated search for this 
signature in 43 months of data from the Fermi Large Area Telescope, building on standard 
methods for gamma-ray line searches. In order to determine the optimal target regions for 
our spectral analysis, we introduced a new adaptive method that takes into account different 
conventional and cuspy DM halo profiles, see Fig. 3. We believe that this method will turn 
out to be useful even in other contexts, essentially whenever spectral features are being looked 
for at the statistical limits of the detector. 

For our toy model, in case of leptonic final states, we find upper limits on the annihilation 
cross-section that are stronger than what can be obtained from dwarf galaxies or collider 
searches, see Figs. 4, 6, 7 and 8. Our limits are still about an order of magnitude too 
weak to constrain annihilation cross-sections naively expected for a thermal relic (assuming 
a standard NFW profile for the DM density and no enhancement of the annihilation rate due 
to substructures). However, future prospects to do so are quite good [26]. 

We also find a weak indication for an internal bremsstrahlung signal, see Fig. 5, corre- 
sponding to a DM mass of ruy^ = 149 ± 4 ~^-^^^ GeV and an annihilation rate of (av) _^ff = 
(6.2 zb 1.5 li'4) X 10~^ cm^ s~^ (we note here that the same signal can also be fitted 
with a conventional gamma-ray line at around 130 GeV and with a cross-section of about 
{av)^^^^f^ ~ 10~^^ cm^ s""*^, and refer to Ref. [81] for a more detailed discussion). After 
taking into account the look-elsewhere effect, the signal significance is about 3.1cr (without, 
it is 4.3(t). We have performed several statistical tests of our method and deem that a purely 
instrumental effect would be a very unlikely explanation for this signal. As stressed above, 
however, such a large radiative annihilation cross-section would be too large to be compatible 
with naive expectations for a thermal relic, at least in our simple scenario and for standard 
cosmological and astrophysical assumptions. In any case, although the observation of a VIB 
or line-like feature from DM annihilation is a fascinating possibility, we caution that more 
data and a much more refined analysis, taking into account all systematics of the LAT, are 
required to reject or confirm this interpretation. 
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A Selection of target regions 

In order to define our target regions, we proceed as follows: We produce a two-dimensional 
equirectangular count map that covers the region defined by ]6] < 90° galactic latitude and 
—90° < i < 90° galactic longitude; the pixel size is A6 = Ai = 1°; each pixel i contains 
the number of gamma-ray events Cj that were measured with energies in the range 1~40 
GeV (note that these Cj are completely unrelated to the Cj introduced in Section 3.3). Since 
gamma-ray fluxes drop rapidly with energy, the countmap produced in this way is completely 
dominated by events with energies close to 1 GeV. We will use this count map as a simple 
but efficient template to estimate the spatial distribution of background events in the sky. 
Since for our spectral analysis we are actually only interested in DM signatures with energies 
larger than 40 GeV, we make the assumption that the fluxes measured at ~ 1 GeV resemble 
the spatial distribution of background events at higher energies reasonably well (instead, one 
could use models for the diffuse emission of the galaxy; this is left for future work). For each 
pixel i, we then calculate the number of expected signal events fii as predicted by (3.1) in the 
energy range 40-300 GeV. This number depends on the adopted dark matter profile (namely 
the value of the inner slope a), and is only defined up to an overall normalization, because we 
leave the annihilation cross-section and the actual annihilation spectrum unspecified at this 
point. Note that the finite angular resolution of Fermi LAT, which is A9 < 0.2° above 40 GeV, 
is neglected. The signal-to-noise level in each pixel can now be estimated as TZi oc ^i/^^/ci 
because a potential signal will likely only be a subdominant perturbation of the background 
fluxes, i.e. fj-i <^ Ci. 

Under the above assumptions, the optimal target region is given by the set of pixels To 
for which the overall signal-to-noise ratio TZ-j-^, defined as 

7^r. = ^m^ , (A.l) 



is maximized. Finding the true To requires in principle a scan over all ~ 2^^*^ possible pixel 
combinations. Since this is unfeasible, we obtain an approximate T by using the following 
simple algorithm: 

1 . We start with an empty set T and include only the one pixel with the largest individual 
signal-to- noise level T^j as seed (this pixel typically lies at the galactic center). 
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Figure 9. Histogram of tlie maximal TS values obtained from a subsampling analysis of the Fermi 
LAT data in the hemisphere pointing towards the galactic anticenter. The blue line shows the theo- 
retical distribution that we used to calculate the look-elsewhere effect. 

2. For each pixel that is not already in T, we calculate how TZj- changes when this pixel 
is added; then all pixels that are found to improve TZj- are added to T at once. 

3. For each pixel in T, we calculate how TZj- changes when this pixel is removed; then all 
pixels for which an increase of 7^7- is found are removed from T at once. 

4. We repeat steps 2 and 3 until the number of pixels in T remains constant. 

5. We remove remaining small regions in T that are not directly connected to the (always 
dominating) region at the galactic center. 

The target region obtained in this way is a very good approximation to the optimal region, 
To — T. Note that the removal of remaining small regions in point 5 does not affect our 
results but merely cleans up the borders of the derived target region, and that the final 
regions are practically independent of the position of the seed pixel in point 1. The regions 
obtained in this way for the different adopted DM profiles are shown in Fig. 3, and will be 
used during our spectral analysis. 

B Details on the statistical analysis 

In order to study the sampling distribution of the test statistic TS (3.3) in absence of a 
signal, we performed a subsampling analysis of Fermi LAT data. To this end, we extracted 
the gamma-ray events measured in the hemisphere pointing towards the anti-galactic center, 
with longitudes \i\ > 90°. Any signal from DM annihilation should be significantly suppressed 
in this direction. From these events we generate 30000 random sample spectra, with the 
Poisson expectation values in each energy bin given by fij = fcj. Here, Cj is the number of 
actually measured events in bin j, and / = 0.13 is adjusted such that the total number of 
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Figure 10. Experimental sensitivity for VIB-like features compared to actual limits. The yellow 
(green) bands show the expected limits at 68% (95%) CL, see text for details. The black solid 
line shows the actually observed limits; these limits are significantly weaker than expected at dark 
matter masses around 150 GeV. The dashed black line shows for comparison the limits obtained when 
removing the data between 115 to 145 GeV from our fits. 



events above 1 GeV in each random sample is ~ 4 x 10^ (in Regl and Reg2 the number of 
events are 5.8 x 10^ and 2.7 x 10^, respectively). In the limit of large event numbers, this is 
equivalent to subsampling the energy distribution of the measured events with replacement. 
In each of these sample spectra, we search for VIB features like discussed above and record 
the largest TS value found. The histogram of the maximal TS values that are obtained in 
this way is shown in Fig. 9. There, we also show the distribution that one obtains when 
selecting the maximum from 4 trials over a x'k=2 distribution. The agreement is very good, 
and we used this distribution when calculating the look-elsewhere effect above. 

In Fig. 10 we show the observed limits (black solid lines) in comparison with the limits 
that are expected at 68% (yellow) and 95% (green) CL. We derived these expected limits 
from 2000 mock data samples that were generated from the null model. In Regl to Reg3, the 
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limits at m^ ~ 150 GeV are significantly weaker than the expectation; this corresponds to 
the large TS values in the left panel of Fig. 5. On the other hand, the relatively strong limits 
at around m^ ^ 100 GeV and m^ ^ 250 GeV are a consequence of the 150 GeV excess, 
which influences the background fits. To illustrate this, we show by the dashed black lines 
the limits that we obtain when removing all data between 115 to 145 GeV (where the excess 
is most pronounced) from the fits; in this case the limits remain in the expected range. 
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